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Abstract. During their main sequence evolution, massive stars can develop convective regions 
very close to their surface. These regions are caused by an opacity peak associated with iron 
ionization. Cantiello et al. (2009) found a possible connection between the presence of sub- 
photospheric convective motions and small scale stochastic velocities in the photosphere of 
early-type stars. This supports a physical mechanism where microturbulence is caused by waves 
that are triggered by subsurface convection zones. They further suggest that clumping in the 
inner parts of the winds of OB stars could be related to subsurface convection, and that the 
convective layers may also be responsible for stochastic excitation of non-radial pulsations. 
Furthermore, magnetic fields produced in the iron convection zone could appear at the surface 
of such massive stars. Therefore subsurface convection could be responsible for the occurrence of 
observable phenomena such as line profile variability and discrete absorption components. These 
phenomena have been observed for decades, but still evade a clear theoretical explanation. Here 
we present preliminary results from 3D MHD simulations of such subsurface convection. 
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1. Introduction 

Hot luminous stars show a variety of phenomena in their photospheres and in their 
winds which still lack a clear physical interpretations at this time. Among these phe- 
nomena are photospheric turbulence, line profile variability (LPVs), discrete absorption 
components (DACs), wind clumping, variable or constant non-thermal X-ray and radio 
emission, chemical composition anomalies, and intrinsic slow rotation. Cantiello et al. 



(2009) argued that a convection zone close to the surface of hot, massive stars, could 
be responsible for various of these phenomena. This convective zone is caused by a peak 
in the opacity due to iron recombination and for this reason is referred as the "iron 
convection zone" (FeCZ). A physical connection may exist between microturbulence in 
hot star atmospheres and a subsurface FeCZ. The strength of the FeCZ is predicted to 
increase with increasing metallicity Z , decreasing effective temperature T and increasing 
luminosity L, and all three predicted trends are reflected in observational data of micro- 



agreement with the results of Cantiello et al. (2009). This suggests that microturbulence 



turbulence obtained in the context of the VLT-FLAMES survey of massive stars (Evans 
|et al.|2005| ). Moreover recent measurements of microturbulence ( Fraser et al.|2010 | are in 



corresponds to a physical motion of the gas in hot star atmospheres. This motion may 
then be connected to wind clumping, since the empirical microturbulent velocities are 
comparable to the local sound speed at the stellar surface. 
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The FeCZ in hot stars might also produce localized surface magnetic fields (Cantiello 



et al. 2009). Such magnetic fields may become buoyant and reach the surface, creating 



magnetic spots. This could explain the occurrence of DACs (discrete absorption compo- 
nents in UV absorption lines) , also in very hot main sequence stars for which pulsational 
instabilities are not predicted. Moreover there may be regions of the upper HR diagram 
for which the presence of the FeCZ influences, or even excites, non-radial stellar pulsa- 
tions. Interestingly stochastic excitation of non-radial pulsations has been recently found 
m massive stars ( |Belkacem et al.||2009| [Degroote et al]|2010[ ). 

The FeCZ could also turn out to directly affect the evolution of hot massive stars. If it 
induces wind clumping, it may alter the stellar wind mass-loss rate. Such a change would 
also influence the angular momentum loss. In addition, magnetic fields produced by the 
iron convection zone could lead to an enhanced rate of angular momentum loss. These 
effects become weaker for lower metallicity, where the FeCZ is less prominent or absent. 

Finally, the consequences of the FeCZ might be strongest in Wolf-Rayet stars. These 
stars are so hot that the iron opacity peak, and therefore the FeCZ, can be directly at 
the stellar surface, or to be more precise at the sonic point of the wind flow (Heger 



&: Langer 1996). This may relate to the very strong clumping found observationally in 
Wolf-Rayet winds ( |Lepine fc Moffat || 1 999| |Marchenko et al.||2006[ ), and may be required 



for an understanding of the very high mass-loss rates of Wolf-Rayet stars (JEichlcr et al. 
1995[|Kato fc Iben||1992[|Heger fc Langer||1996). 



2. Simulations of subsurface convection 

Convection in the FeCZ is relatively inefficient: the transport of energy is dominated 
by radiation, which usually accounts for more than 95% of the total flux. This region of 
the star is very close to the photosphere, above which strong winds are accelerated. The 
continuous loss of mass from the stellar surface moves the convection region to deeper 
layers, revealing to surface material that has been processed in the FeCZ. In rotating 
stars, the associated angular momentum loss might also drive strong differential rotation 
in the region of interest. It is clear that, under these circumstances, the mixing length 
theory can only give a qualitative picture of the convective properties in these layers. In 
order to study the effects induced by the presence of subsurface convection at the stellar 
surface in a more quantitative way, we perform local 3D MHD calculations of convection. 
In these simulations we can vary the relative importance of the background radiative flux 
and include the effects of rotation and shear, in order to model conditions as similar as 
computationally permitted to subsurface convection in OB stars. 



2.1. Computational model 
The setup is similar to that used by Kapyla, et al. ( 2008 ) . A rectangular portion of a star is 



modeled by a box situated at colatitude 6. The dimensions of the computational domain 
are {L x , L yi L z ) = (5, 5, 5)d where d is the depth of the convectively unstable layer. Our 
(x,y,z) correspond to (0, 0, r) in spherical polar coordinates. The box is divided into 
three layers, an upper cooling layer, a convectively unstable layer, and a stable overshoot 
layer (see below) . The following set of equations for compressible magnetohydrodynamics 
is being solved: 



dA 

~dt 



Sx 



dA 

dy 



= U x B rj^oJ - SAyX, 



(2.1) 



(2.2) 
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where D/Dt = d/dt + (U + Uo) ■ V, and Uo — (0, 5a;, 0) is an imposed large-scale shear 
flow in the y-direction. The magnetic field is written in terms of the magnetic vector 
potential, A, with B = V x A, J = /< ] V x B is the current density, po is the vacuum 
permeability, 77 and v are the magnetic diffusivity and kinematic viscosity, respectively, 
K is the heat conductivity, p is the density, U is the velocity, g = —gz is the gravitational 
acceleration, and i~2 = fto(— sm8, 0, cos 9) is the rotation vector. The fluid obeys an ideal 
gas law p — (7 — l)pe, where p and e are pressure and internal energy, respectively, 
and 7 = cp/cy = 5/3 is the ratio of specific heats at constant pressure and volume, 
respectively. The specific internal energy per unit mass is related to the temperature via 
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Figure 1. Simulation of subsurface convection. Starting from the top, the box is divided into 
three layers: a radiative layer with an upper cooling boundary, a convectively unstable layer and 
another stable layer at the bottom. The snapshot shows values of vertical velocity. The plane 
below the box shows the vertical velocity field at the lower boundary of the convective layer. 
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Figure 2. Left: evolution of the magnetic field as function of time (in code units) in three 
different runs with different physics. The value on the y-axis is the logarithm of the ratio between 
the magnetic field and the equipartition value (the value corresponding to equipartition between 
kinetic and magnetic energy). Dynamo action reaching equipartition is found in simulations 
including both rotation and shear. Right: energy spectrum of the velocity field in our run of 
subsurface convection. On the x-axis is the spatial wavenumber k, on the j/-axis the z coordinate 
of the computational domain in gridpoints. The color bar on the right indicates the value of the 
energy per wavenumber E(k). The location of the convection zone in the z-coordinate is shown 
by the dashed lines, the rest of the domain is radiative. Just above and below the convection 
zone kinetic energy is associated with overshooting, while at further distances the transport of 
energy is facilitated by gravity waves. 



e = c\/T . The rate of strain tensor S is given by 

Sy = U U hi + U i,i) - ¥v V ■ U ( 2 - 5 ) 



The last term of Eq. (2.4) describes cooling at the top of the domain. Here t(z) is a 
cooling time which has a profile smoothly connecting the upper cooling layer and the 
convectively unstable layer below, where r —¥ oo. 

The positions of the bottom of the box, bottom and top of the convectively unstable 
layer, and the top of the box, respectively, are given by (z 1; z 2 , z 3 , z±) = (—0.5, 0, 1, 4.5)<i, 
where d is the depth of the convectively unstable layer. Initially the stratification is 
piecewise polytropic with polytropic indices (7711,7712,7713) = (3,0.9,3), which leads to 
a convectively unstable layer between two stable layers. The cooling term leads to a 
stably stratified isothermal layer at the top. All simulations with rotation use 8 = 0° 
corresponding to the north pole. 

Stress-free boundary conditions are used in the vertical (z) direction for the velocity, 

U x<z = U y . z =U Z = 0, (2.6) 

where commas denote partial derivatives, and perfect conductor conditions are used for 
the magnetic field, i.e. 

B x , z = B y>z = B Z =0. (2.7) 

In the x and y directions periodic boundary conditions are used. The simulations were 
performed with the Pencil CODE[f] which uses sixth order explicit finite differences in 
space and third order accurate time stepping method. Resolutions of up to 128 2 x 256 
mesh points were used. 



t http : //pencil-code . googlecode . com/ 
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2.2. Preliminary results 

As a preliminary study we performed low resolution simulations (128x128x256), where 
the density contrast between the bottom of the convective layer and the top of the domain 
is only ^20. This is about ten times smaller than in the case of the FeCZ. Moreover the 
ratio of the convective to radiative flux is about 0.3, higher than in the FeCZ case. 
This is because smaller values of convective flux result in steady convection at the low 
resolution of these preliminary runs. Therefore, at this stage, the velocities of convective 
motions cannot directly be compared to the velocities obtained by mixing length theory. 
However already in these preliminary runs we could follow the excitation and propagation 
of gravity waves above the convective region. In the right panel of Fig. [2] we show the 
kinetic energy spectrum (as function of the spatial wavenumber k) in the horizontal 
plane, as function of depth. The maximum of energy is found in the convective region for 
those wavenumbers k corresponding to the number of resolved convective cells (about 5 
along one of the horizontal directions). Energy is also transported up to the top layer by 
gravity waves, where the maximum of the energy is deposited in those wavelengths that 
are resonant with the scale of convective motions, as predicted, for example, by |Goldreich| 



& Kumar ( 1990 



Kapyla et al. ( 2008 1 found excitation of a large scale dynamo in simulations of turbulent 



convection including rotation and shear. Our computational setup is very similar, so it's 
not surprising that we can confirm this result. Dynamo action reaching equipartition is 
found in our simulations that include shear and rotation (see left panel of Fig. [2| , with 
magnetic fields on scales larger than the scale of convection. 

2.3. Discussion 



The connection found by Cantiello et al. ( 2009 1 between the presence of sub-photosphcric 
convective motions and microturbulence in early-type stars is intriguing. The recent 



identification of solar-like oscillations in hot massive stars (Belkacem et al. 2009 Degroote 
et al.pOlO I and further measurements of microturbulence ( Fraser et al.||2010" ) also point 



toward a picture in which the FeCZ influences surface properties of OB stars. 

We performed 3D MHD simulations of convection to investigate the excitation and 
propagation of gravity waves above a subsurface convection zone. Analytical predictions 



of Goldreich & Kumar ( 1990 ) on the spatial scale at which the maximum of energy is 



injected in gravity waves seem to be confirmed in our preliminary calculations. Further 
investigation is required to understand if the subsurface convection expected in OB stars 
excites gravity waves of the required amplitude to explain the observed microturbulence 
in massive stars. In particular we need higher resolution to increase the Reynolds number 
of our simulations and be able to decrease the ratio of convective to radiative flux, 
which appear to be an important parameter in determining the convective velocities 
( Brandenburg et al.||2005[ ). 

Magnetic fields reaching equipartition values are found in simulations of turbulent 
convection if rotation and shear are present (Kapyla et al. 2008). Since massive stars are 
usually fast rotators, it could be that the interplay between convection, rotation and shear 
is able to drive a dynamo in OB stars. Indeed our simulations of subsurface convection 
including rotation and shear show the excitation of dynamo action, with magnetic fields 
reaching equipartition. This means that fields up to a kG could be present in the FeCZ. 
These magnetic fields might experience buoyant rise and reach the surface of OB stars, 
where they could have important observational consequences. In particular it has already 
been suggested that the discrete absorption components observed in UV lines of massive 
stars could be produced by low amplitude, small scale magnetic fields at the stellar 



surface (Kaper & Henrichs 1994). Further study is needed to investigate the amplitude 
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and geometry of magnetic fields reaching the stellar surface, as well as their effects on 
the photosphere. 
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